home *** CD-ROM | disk | FTP | other *** search
/ Ahoy 1989 February / Ahoy_Magazine_89-02_1989_Double_L.d64 / Dio Solver (.txt) < prev    next >
Commodore BASIC  |  2022-10-26  |  4KB  |  115 lines
  1. 1 rem =================================
  2. 2 rem
  3. 3 rem       diophantine solver
  4. 4 rem        rupert report #62
  5. 5 rem ---------------------------------
  6. 6 rem  solve equations of the form
  7. 7 rem   ax + by = c  for x and y
  8. 8 rem  where all terms are integers
  9. 9 rem =================================
  10. 10 gosub 100   :rem initialize
  11. 20 gosub 300   :rem get a,b,c
  12. 30 gosub 1000  :rem make a & b positive
  13. 40 gosub 2000       :rem convert a/b to  continued fraction and get gcf(a,b)
  14. 50 gosub 3000 :rem table of convergents
  15. 60 gosub 4000 :rem find t for x,y > 0
  16. 70 gosub 5000 :rem find x,y for given t
  17. 90 end
  18. 99 rem --------------------------------
  19. 100 rem -- initialize & input variables
  20. 110 dim p(20),q(20),f(20)
  21. 120 def fndiv(b)=int(a/b)   :rem integer division
  22. 130 def fnmod(b)=int(a-b*fndiv(b))  :rem mod function
  23. 140 def fnii(z)=(z=int(z))  :rem true if z is integer
  24. 150 def fnmi(x)=sgn(x)*int(abs(x))  :rem make integer
  25. 160 def fnx(t)=c*q0-snb*b*t    :rem x(t)
  26. 170 def fny(t)=a*t-snb*c*p0    :rem y(t)
  27. 175 print "[147]"
  28. 180 print: print "enter non-0 integer coefficients a,b,c"
  29. 190 print "for equation  ax + by = c:"
  30. 200 input a,b,c
  31. 210 if not (fnii(a) and fnii(b) and fnii(c)) then print"integers only": goto 180
  32. 220 if a=0 or b=0 or c=0 then print"non-zero integers only": goto 180
  33. 230 print "[147]equation is ";a;"* x + ";b;"* y = ";c
  34. 240 print
  35. 250 a0=a: b0=b: c0=c: rem save originals
  36. 260 return
  37. 300 rem -- convert to positive integers
  38. 310 snb=1   :rem sign of b
  39. 320 if a<0 then a=-a: b=-b: c=-c    :rem make a>0
  40. 330 if b<0 then b=-b: snb=-1   :rem make b>0
  41. 340 a1=a: b1=b: c1=c   :rem save working values
  42. 350 return
  43. 360 rem --------------------------------
  44. 1000 rem   convert a/b to continued fraction f(n) and find gcf
  45. 1010 rem  a/b = f(n) remainder m
  46. 1020 n=1  :rem number of terms
  47. 1030 f(n)=fndiv(b) :rem integer division
  48. 1040 m=fnmod(b)    :rem remainder
  49. 1050 if m=0 then 1080  :rem done
  50. 1060 a=b: b=m: n=n+1
  51. 1070 goto 1030
  52. 1080 if fnii(n/2) then 1110  :rem n/2 is integer::n is even
  53. 1090 f(n)=f(n)-1    :rem convert to even number of terms
  54. 1100 f(n+1)=1: n=n+1
  55. 1110 gcf=b  :rem final divisor is gcf of a,b
  56. 1120 a=a1: b=b1: c=c1 :rem restore working values
  57. 1130 rem -----next line prints continued fraction for a/b
  58. 1140 print a; "/"; b; "= [";: for k=1 to n: print f(k);: next k: print "]"
  59. 1150 print "  gcf ("; a; ","; b; ") ="; gcf: print
  60. 1160 return
  61. 1170 rem -------------------------------
  62. 2000 rem --reduce eqn by gcf if possible
  63. 2010 if gcf=1 then 2050
  64. 2020 if fnii(c/gcf) then goto 2040
  65. 2030 print "no solutions.  c is not divisible by gcf of a and b": end
  66. 2040 a=a1/gcf: b=b1/gcf: c=c1/gcf
  67. 2050 return
  68. 2060 rem -------------------------------
  69. 3000 rem -- table of convergents of a/b
  70. 3010 p(0)=0: p(1)=1
  71. 3020 q(0)=1: q(1)=0
  72. 3030 rem print "convergents p,q:"
  73. 3040 for k=2 to n+1
  74. 3050 p(k)=f(k-1)*p(k-1)+p(k-2)
  75. 3060 q(k)=f(k-1)*q(k-1)+q(k-2)
  76. 3070 rem print p(k),q(k)
  77. 3080 next k
  78. 3090 p0=p(n): q0=q(n)
  79. 3100 return
  80. 3110 rem ------------------------------
  81. 4000 rem -- show results
  82. 4010 print: print "x and y as functions of t :"
  83. 4020 print " (t = 0, +/-1, +/-2, ...)":  print
  84. 4030 sn$=") - (": if snb=-1 then sn$=") + ("
  85. 4040 print "   x = ("; c*q0; sn$; b; "* t )"
  86. 4050 print "   y = ("; a; "* t "; sn$; c*p0; ")"
  87. 4060 print: print "values of t for x & y > 0:"
  88. 4070 tx=snb*c*q0/b
  89. 4080 tx$=" < ": if snb=-1 then tx$=" > "
  90. 4090 ty=snb*c*p0/a: ty$=" > "
  91. 4100 print " for x>0: t"; tx$; tx
  92. 4110 print " for y>0: t"; ty$; ty
  93. 4120 if snb=-1 then 4200
  94. 4130 rem  find t for t<tx and t>ty
  95. 4140 if tx < ty then print "no integer solutions for x>0 and y>0": goto 4230
  96. 4150 ta=int(tx): if ta=tx then ta=ta-1
  97. 4160 tb=int(ty)+1
  98. 4170 if ta<tb then print "no integer solutions": goto 4230
  99. 4180 if ta=tb then print "only solution is t = "; ta: goto 4230
  100. 4190 print "t is any integer from "; tb; " to "; ta: goto 4230
  101. 4200 rem  find t for t>tx and t>ty
  102. 4210 mn=int(tx)+1: if ty>tx then mn=int(ty)+1
  103. 4220 print "t is any integer greater than or equal to "; mn
  104. 4230 print
  105. 4240 return
  106. 4250 rem ------------------------------
  107. 5000 rem --find x,y for given t
  108. 5010 input "what value of t";t
  109. 5020 print " x is ";fnx(t)
  110. 5030 print " y is ";fny(t)
  111. 5040 print: input"more [y,n]y[157][157][157]"; a$
  112. 5050 if a$="n" then return
  113. 5060 print "[145]              ": print"[145]";
  114. 5070 goto 5010
  115.